/*! 
 ***************************************************************************
 * \file BuckleyLeverett1D.cpp
 * Two phase flow in porous media in 1D.
 *
 * \par Description
 * This code solves the two phase flow in porous media in 1D, the
 * equations are posed in pressure-saturation formulation:
 * 
 * These equation are solved using the IMPES algorithm.
 * The domain of study is the interval \f$ x \in [0,300 m] \f$, and is
 * initially filled with oil. Then, water is injected on the left side
 * and we study the traveling front between oil and water.
 *
 *  \par Compiling and running
 *  Modify the variables BASE and BLITZ in the file \c tuna-cfd-rules.in 
 *  according to your installation and then type the next commands:
 *  \verbatim 
   % make
   % ./BuckleyLeverett1D  \endverbatim
 ***************************************************************************
 *  \author    Luis M. de la Cruz [ Mon Mar 17 08:47:37 GMT 2010 ]
 ***************************************************************************
 */

#include <iostream>
#include <cmath>
#include "Meshes/Uniform.hpp"
#include "Storage/DiagonalMatrix.hpp" 
#include "Equations/TwoPhaseEquation.hpp"
#include "./BL/BLIP1.hpp"
#include "./BL/BLES1.hpp"
#include "Solvers/TDMA.hpp"
#include "Solvers/explicit.hpp"
#include "Utils/inout.hpp" 
#include "Utils/num_utils.hpp"
#include "Utils/GNUplot.hpp"
using namespace Tuna;
  
// Scalar field in one dimension
typedef TunaArray<double,1>::huge ScalarField1D;

int main() 
{
    double length, dt, Tmax, tolerance; 
    int num_nodes, frequency;
    double porosity, permeability, mu_w, mu_o, Srw, Sro;
    double injection, pres_out, Sw_in;

// Reading data from the input file
    std::ifstream input_cin ("inputBL1D"); 
    input_cin >> length
	      >> num_nodes
	      >> porosity
	      >> permeability
	      >> mu_w
	      >> mu_o
	      >> Srw
	      >> Sro
	      >> injection
	      >> pres_out
	      >> Sw_in
	      >> dt
	      >> Tmax
	      >> frequency;
    input_cin.close();

    double Gamma = permeability;

    StructuredMesh<Uniform<double, 1> > mesh(length, num_nodes);
    float dx = mesh.getDelta(X);
    mesh.print();

    ScalarField1D p ( mesh.getExtentVolumes() );
    ScalarField1D Sw ( mesh.getExtentVolumes() ); 
    ScalarField1D p_n( mesh.getExtentNodes()   );
    ScalarField1D Sw_n( mesh.getExtentNodes()   );

    // Initial conditions
    Sw = 0;
    p = pres_out; 
 
    // Boundary condition
    Sw(0) = Sw_in;  

    NumUtils::interpolateToNodes(p_n, p);
    InOut::writeToFile(p_n, 0, "./DataBL1D/pres.", dx);
    NumUtils::interpolateToNodes(Sw_n, Sw);
    InOut::writeToFile(Sw_n, 0, "./DataBL1D/satu.", dx);

        
#ifdef WITH_GNUPLOT
    GNUplot plotter;
    int pause; 
    plotter("set grid");
    /*
    plotter("p \"./DataBL1D/pres.00000\" title \"Pressure\"w lp");
    std::cout << "\n >---> Press <1> and then <enter> to continue \n";
    std::cin >> pause;
    plotter("p \"./DataBL1D/satu.00000\" title \"Saturation\"w lp");
    std::cout << "\n >---> Press <1> and then <enter> to continue \n";
    std::cin >> pause;
    */
#endif

    SparseMatrix< Diagonal< double, 1> > A(num_nodes); 
    //    DiagonalMatrix< double, 1> A(num_nodes); 
    ScalarField1D              b(num_nodes);

    TwoPhaseEquation< BLIP1<double, 1> > pressure(p, A, b, mesh.getDeltas());
    pressure.setDeltaTime(dt); 
    pressure.setPermeability(permeability);
    pressure.setPorosity(porosity);
    pressure.setSrw(Srw);
    pressure.setSro(Sro);
    pressure.setViscosity_w(mu_w);
    pressure.setViscosity_o(mu_o);
    pressure.setInjection(injection);
    pressure.setNeumann  (LEFT_WALL, -injection * mu_w / permeability);
    pressure.setDirichlet(RIGHT_WALL);
    pressure.setSaturation(Sw);
    pressure.print();

    //
    // In these case, the boundary conditions are calculated in
    // the numerical scheme (the adaptor BLES1).
    //
    TwoPhaseEquation< BLES1<double, 1> > saturation(Sw, A, b, mesh.getDeltas());
    saturation.setDeltaTime(dt); 
    saturation.setPermeability(permeability);
    saturation.setPorosity(porosity);
    saturation.setSrw(Srw);
    saturation.setSro(Sro);
    saturation.setViscosity_w(mu_w);
    saturation.setViscosity_o(mu_o);
    saturation.setInjection(injection);   
    saturation.applyBounds(1, Sw.ubound(firstDim)-1);// just to move the indexes
    saturation.setPressure(p);
    saturation.print();

    ofstream error_fp("./DataBL1D/error"), residual_fp("./DataBL1D/residual");
    int iteration = 1, idprint = 1;
    double t = dt;
    double error = 1.0;
    double residual = 1.0;
    int pausa;

    int bi = saturation.get_bi(), ei = saturation.get_ei();

    frequency *= (86400 / dt);

    while (t <= Tmax) {

      pressure.calcCoefficients();
      Solver::TDMA1D(pressure);
      residual = pressure.calcResidual();
      error = pressure.calcErrorL1();
      error_fp << iteration << "\t" << error << std::endl;
      residual_fp << iteration << "\t" << residual << std::endl;	
      pressure.update();

      saturation.calcCoefficients();
      Solver::solExplicit(saturation);
      error = saturation.calcErrorL1();
      saturation.update();

      /* *
      // Flux limiter
      for(int i = 0; i <= (ei+1); ++i) {
	if ( Sw(i) > Sw_in) Sw(i) = Sw_in;
	if ( Sw(i) < Srw  ) Sw(i) = Srw;
      }
      saturation.updatePhi(Sw);
      /* */

      if ( !(iteration % frequency) ) {
	NumUtils::interpolateToNodes(p_n, p);
	InOut::writeToFile(p_n, 1, "./DataBL1D/pres.",dx);

	NumUtils::interpolateToNodes(Sw_n, Sw);
	InOut::writeToFile(Sw_n, 1, "./DataBL1D/satu.",dx);

	/* *
#ifdef WITH_GNUPLOT
 	plotter("p \"./DataBL1D/satu.00001\" title \"Saturation\" w lp");
 	std::cin >> pause;
#endif
	/* */

	std::cout << "\n ---- ID print => " << idprint
		  << "\t - Iteration => " << iteration 
		  << "\t Step = " << t
		  << "\n ---- Res = " << residual
		  << "\t Err = " << error;
	std::cout.flush();
	++idprint;
      }

      ++iteration;
      t += dt;

    }

#ifdef WITH_GNUPLOT
	plotter("p \"./DataBL1D/pres.00001\" title \"Pressure\" w lp");
	std::cout << "\n \n Press 1 and enter to see the saturation ";
	std::cin >> pause;
	plotter("p [][0:1] \"./DataBL1D/satu.00001\" title \"Saturation\" w lp");
	std::cout << "\n \n Press 1 and enter to finish ";
	std::cin >> pause;
#endif	
    /** **/


    return 0;
}




